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Abstract. A peculiarity of tlie iiydrodynamic Navier- Stokes equations for a granular gas is the modification of the Fourier 
law, with the presence of an additional contribution to the heat flux that is proportional to the density gradient. Consequently, 
the constitutive relation involves, in the case of a one-component granular gas, two transport coefficients: the usual (thermal) 
heat conductivity and a diffusive heat conductivity. A very simple physical interpretation of this effect, in terms of the mean 
free path and the mean free time is provided. It leads to the modified Fourier law with an expression for the diffusive 
Fourier coefficient that differs in a factor of the order of unity from the expression obtained by means of the inelastic 
Boltzmann equation. Also, some aspects of the Chapman-Enskog computation of the new transport coefficients as well as 
of the comparison between simulation results and theory are discussed. 
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INTRODUCTION 

Granular fluids are of increasing interest to several scientific communities. This is partially due to their great relevance 
in many industrial issues of agricultural, pharmaceutical and chemical significance for packing and transport of grains. 
Moreover, they raise a number of important conceptual challenges. A class of these questions, widely investigated 
in the last two decades or so, involve the form and conditions for a hydrodynamic description, similar to the one for 
normal, molecular fluids [1, 2, 3, 4]. In most of the cases, the simplest model of a granular gas at the particle level 
of description, an ensemble of inelastic hard spheres or disks [5], is considered. This model has proven to be able to 
describe at a qualitative level many of the peculiarities exhibited by real granular gases. 

The Navier-Stokes equation for the energy of a granular gas exhibits two main peculiarities when compared with 
that for normal fluids. First, there is a source term associated to the energy dissipation in collisions. Its origin can be 
easily understood on the basis of mean field arguments. The other differential feature is the modification of the Fourier 
law for the heat flux, including a new term that couples it to the density gradient. This new contribution to the heat flux 
was predicted by kinetic theory methods [6, 7, 8], and later on measured in molecular dynamic simulations [9]. Also, 
some implications of its existence have been observed in experiments [10, 11, 12, 13]. The theoretical formal reasons 
why this density gradient contribution to the heat flux is present in granular fluids and absent in normal ones has been 
investigated. They are related to the time reversal invariance of both the dynamics of the system and the equilibrium 
distribution function of normal fluids [14]. 

On the other hand, although some physical interpretations of the modification of the Fourier law in granular gases 
have been given, it seems instructive to provide simple intuitive arguments leading to the new term in the Fourier law. 
In particular, one of the questions addressed here is to show that mean free path and collision time arguments are 
enough as to predict the presence of the density gradient contribution to the heat flux for granular fluids, in a similar 
way as they predict all the transport coefficients in molecular gases, aside from a numerical factor 

Another aim of this paper is to discuss a new approximation method to compute Navier-Stokes transport coefficients 
of dilute gases starting from their Green-Kubo representation. The method is both simple and transparent. It will be 
illustrated here for the two coefficients appearing in the modified Fourier law. In the elastic limit the results are the 
same as those derived by the usual Sonine expansion to lowest order, while in the inelastic limit both results are very 
close. 

The remaining of the paper is organized as follows. In the next section, the derivation of the generalized Fourier 
law using linear response theory in the context of the Boltzmann equation is shortly reviewed. The relevance of the 
hydrodynamic part of the spectrum of the linearized Boltzmann operator is emphasized. The formal expressions for the 
transport coefficients are explicitly evaluated in Sec. by using an approximation related with the non-hydrodynamic 
part of the spectrum of the linear Boltzmann operator The results are shown to be practically equivalent to those that 
have been obtained by expanding in Sonine polynomials keeping only the lowest order 



In Sec. a much simpler derivation of the modified Fourier law is carried out. It is based on the concepts of mean free 
time and mean collision time in a dilute gas, generalizing the arguments used in elementary kinetic theory of molecular 
gases [16]. The origin of both contributions to the heat flux, as well as the way in which the inelasticity modifies the 
usual term proportional to the temperature gradient, show up in a transparent way. Finally, the last section contains 
a short discussion of some macroscopic effects associated with the new term in the expression of the heat flux. The 
specific state considered corresponds to a series of experiements carried out by different groups. 



LINEARIZED BOLTZMANN EQUATION AND HYDRODYNAMICS 



The system considered is a dilute granular gas of smooth inelastic hard spheres (d = 2) or disks (d = 3) of mass m 
and diameter a. The coefficient of normal restitution a will be taken as a velocity-independent constant, although the 
theory has also been extended to more realistic models in which it depends on the relative velocity of the colliding 
particles [15]. The one-particle distribution function of the system /(r,v,f) obeys the inelastic nonlinear Boltzmann 
equation [17, 18]. This equation has a particular solution describing the homogeneous cooling state (HCS), 



fH{v,t)=nv/{t)x{c), vo; 



(1) 



where n is the uniform number density, X is an isotropic function of the scaled velocity c, and Th is the uniform 
temperature of the system. It decays monotonically in time according with the Haff law, 



(2) 



Approximated expressions for the distribution x{c) and for the cooling rate i^// have been obtained by expanding the 
function in Sonine polynomials and keeping only the lowest order polynomial [18, 19]. The approximated distribution 
function has the form 
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The expression for the cooUng rate is: 
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Here £ = («a^ * ) 'is proportional to the mean free path of the gas. Suppose now a smaU perturbation around the 
HCS, and define 5f by 



/(r,v,0=/ff(v,0 + 5/(r,v,0, |5/(r,v,0 « |/h(v,0|. 



(7) 



To eliminate the time dependence associated to the HCS, it is convenient to introduce dimensionless length / and time 
scales s defined by 
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respectively. The time scale s is proportional to the accumulated average number of colUsions per particle. The 
dimensionless form of the deviation of the distribution function is 



5x{l,c,s) = n ^vi{t)5f{r,v,t). 



(9) 



This function obeys the linear equation [20, 21 J 



^ + c--^\8x{l,c,s)=A{c)x{l,c,s). (10) 



,ds dl 
The linear Boltzmann operator A is given by 

A(ci ) = y ^;C2 To (ci , C2) ( 1 + /'i2);t (C2) - y ^ • Cl . 

where To{c\ , C2) is the dimensionless binary collision operator for inelastic hard spheres or disks, 

To{cuC2) = J dd@{ci2-B)cu-d[a-^ba\cuC2)-l]. 



(11) 



(12) 



Here Cu = C\—C2, da is the solid angle element for the unit vector a, & is the Heaviside step function, and b^^ {ci , C2) 
is an operator changing all the functions of Ci and C2 to its right by the same functions of the precoUisional velocities 
c\ and C2, given by 



C*2 = b„^C2=C2+^^{a-Cu)0. (13) 

Solutions to Eq. (10) are sought in a Hilbert space defined by the scalar product 

{g\h) = jdcx-\c)g*{c)h{c), (14) 

with g* [c) being the complex conjugate of g{c). It is then useful to consider the homogeneous eigenvalue problem 

A(c)^p(c)=A/3^p(c). (15) 

The solutions of this equation corresponding to the infinite wavelength limit of the hydrodynamic equations are given 
by [20, 21] 

Ai=0, A2 = |, ^3 = -y, (16) 

^Ac)=X{c) + ~[cx{c)], §2^ = -^, ^3(c) = -^-[c;c(c)]. (17) 

The eigenvalue X2 is d-fo\d degenerated. The operator A(c) is not Hermitian and the eigenfunctions are not 
orthogonal. This leads to introduce a set of functions 1^, that be orthogonal to the above eigenfunctions. A convenient 
choice is 

|l(c)=x(c), Uc)=CX{c), 13(e) =(^^ + l^;^(c); (f^|<^/3,)=5/3,^,. (18) 

The general solution of the linearized Boltzmann equation in the Hilbert space can be formally written in the Fourier 
representation as 

8x = ZH{kM?{1^,c) + 5x"'{k,c,s), (19) 
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where the sum extends over the d + 2 hydrodynamic modes of the operator A—ikc and dx'" contains all the other 
"microscopic" modes. The hydrodynamic modes for finite k of A — ik- c are defined as those eigenvalues that are 
continuously connected as functions as k to those given in Eqs. (16). 

It is now assumed that the hydrodynamic part of the spectrum dominates for long times and small gradients. Then, 
the microscopic part can be neglected in Eq. (19). Then, the formal solution of Eq. (10) to first order in k can be 
expressed as [20] 



5x{k,c,s)=e^'^5x{k,c,Q)-Y^ r ds' e''i^-h)c^p{c) 
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■ikaii{k,s). (20) 



Moreover, the coefficients ap in the above expression can be identified as 



ap{k,s):^{^p\5x{k,s)). (21) 

These coefficients can be expressed as functions of the hydrodynamic fields density, velocity, and temperature. In this 
way, a formal expression for dx valid to first order in the gradients (k) is obtained and, using it, the Navier-Stokes 
expressions for the heat flux and the pressure tensor. They are a generalization of the well-known Green-Kubo formulas 
for molecular gases. The details of the calculations are given in ref. [20]. In the next section the expression for the heat 
flux flux will analyzed in some detail. 



GREEN-KUBO EXPRESSION FOR THE HEAT FLUX 

Following the procedure sketched in the previous section, it is obtained that the heat flux q{r,t) to Navier-Stokes order 
has the form 

^ = -(cVT-^iVn, (22) 

where K is the (thermal) heat conductivity and ji a transport coefficient vanishing in the elastic Umit and that is referred 
to as the diffusive heat conductivity. Their formal expressions are 

K = nm£vo{t)K{s), n =m£vl{t)Ji{s), (23) 

Kis)^^Jdci:ic)-^3ic,s), (24) 

Ji{s) = ^Jdci:{c)-^i{c,s), (25) 

with 

E(c)^(c^-^)c. (26) 

«I>i(c,*)= rrf//^W^i(c)c + 24.3(c,*), (27) 
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(28) 



Here (c) and ^2 (c) are the eigenfunctions given by Eqs. (17). As usual, k{s) and ]1 (s) are expected to reach steady 
plateau values for large enough s, when the hydrodynamic description is accurate. Both transport coefficients have 
been evaluated in the first Sonine approximation [7, 22]. Here an alternative approximation will be discussed [23]. It 
consists in treating T.{c)x{c) as an eigenfunction of the operator adjoint of A, A+, i.e. it is considered that 

A+{c)^c)x{c) :^ l5^{c)x{c), (29) 

and A5 is obtained by multiplying this equation by Cx and afterwards integrating over the velocity c. The results reads 

J 4/(a) , Co(«) , 3Co(«) 

{d + 2)a2{a)^a2{a)^ 2 ' ^^^^ 

7(a) = ^ ^ {l6{2 + d)il-a) + a2{a)[70 + 47d-3 (34 + 5d)a]}. (31) 

32v2c/r {d/2) 

Note that assuming Eq. (29) leads to an exponential decay of the time correlation function involved in the calculation 
of K (see Eqs. (24) and (28)). With this approximation, the long time limit value of the thermal heat conductivity k{s) 
is given by 



(32) 




FIGURE 1. Diffusive heat conductivity of a dilute granular gas of hard disks (left side) and spheres (right side) as a function of 
the coefficient of normal restitution a. The solid lines are the results found in the first Sonine approximation, while the dashed lines 
have been obtained using Eq. (34). 



Evaluation of the velocity integral using Eqs. (17) and (3) yields 

~,Jl±m±2a^_ (33) 
2(2|A5|-Co) 

By proceeding in an analogous way, the expression for the diffusive heat conductivity can be evaluated giving 

~^^2~.-^^±^m^. (34) 

4|A5| 

As mentioned above, these coefficients have also been computed in the first Sonine approximation [7, 22]. The results 
obtained in that way are almost indistinguishable from Eqs. (33) and (34), as it can be seen in Fig. 1 for the coefficient of 
diffusive heat conductivity /x. Similar results are obtained for the coefficient of thermal heat conductivity K. It is worth 
to mention that a modified Sonine expansion, in which the Gaussian is replaced by the (approximated) distribution 
of the HCS, has also been employed to compute the transport coefficients of a dilute granular gas [24]. The results 
reported are not equivalent to Eqs. (33) and (34), although they are very close for a > 0.65. 



ELEIVIENTARY DERIVATION OF THE IVIODIFIED FOURIER LAW 

Suppose a three-dimensional granular gas in an arbitrary time-dependent state, but with macroscopic gradients only in 
the z-direction. It has already been mentioned that there is a homogeneous reference state, the homogeneous cooling 
state, in which the temperature of the system decays monotonically in time according to the Haff law. The interest here 
is to derive an expression for the heat flux occurring in the system in the direction of the gradients when they are small. 
Since the heat flux is by definition the flux of internal energy that is not associated to a macroscopic flux of mass, it 
will be assumed that the average flux of particles vanishes everywhere in the system. 

Consider a plane z = constant inside the gas (see Fig. 2). In the dilute Umit, the flux of particles crossing it in the 
direction of increasing z can be estimated as 

^^+\z,t)^\n^-\z,M-\z,t). (35) 

Here n!^^\z,t) and v^^' (z, f ) are the number density and average velocity, respectively, of the granular gas, just below 
the considered plane. Upon writing the above expression, it has been assumed that the velocity distribution in that 



FIGURE 2. Sketch of the way in which the fluxes through a plane z = constant in the interior of the granular gas are computed. 



region can be accurately approximated by an isotropic function. Similarly, the flux of particles through the same plane 
but in the decreasing z-direction is 

« inW(z,OvW(z,j). (36) 

where now (z, f ) and v^"*") {z-t) refer to values just above the plane z = constant. 
The requirement that there is no net flux of particles imphes that 

^(+)(z,0 = ^(-)(z,0 = ^{z,t) « \n{zmz,t). (37) 

Next, the average value, e(+'(z,f) of the energy being carried out by a particle crossing the z-plane at time t in the 
positive direction will be estimated. Typically, the last collision suffered by the particle before crossing the plane was 
at a distance of the order of the local mean free path of the gas, X{z,t). Moreover, that collision took place at a time 
t—x, where x is the mean free time, i.e. the average time between consecutive collisions of a given particle. In the 
spirit of mean free path theories, it is assumed that the effect of the collision can be modeUed as the thermaUzation 
of the velocity distribution of the particles to the local temperature. Consequently, the average energy of a particle 
crossing the plane in the direction of increasing z is approximated by 

e^+\z,t)^^-T{z-X,t-x). (38) 

In the same way, the average energy carried out by a particle crossing the plane in the direction of decreasing z is 
estimated as 

e^-\z..t)^^-T{z + X,t-x). (39) 
Consequently, the neat heat flux through the z-plane is given by 

(lz{z,t)^^n{z, t)v{z, t)[T{z-?i,t-x)-T{z + ?i,t-x)]. (40) 

The aim now is to analyze the above expression in the limit of small gradients of the hydrodynamic fields, namely 
the number density and the temperature. This requires some care when approximating the temperatures difference 
appearing on the right hand side of the equation. It is convenient to expand the temperature difference as 

T{z-X,t-x)-T{z + k,t-T)p^T{z-k,t)-T{z + k,t)-x^[T{z-X,t)-T{z + k,t)]. (41) 



This approximation requires that the temperature profile of the system changes very little over distances of the order of 
the mean free path X for times of the order of the collision time T. Given that the difference r(z — A , f ) — T (z + A , f ) is 
at least of first order in the gradients, its time derivative must be computed without introducing any additional gradient 
operator. Therefore, 

jJ{zT^,T) « -r(zTA,0C(2TA,r), (42) 

where l^{z,t) is the cooling rate to zeroth order in the gradients, i.e. the cooling rate of the homogeneous cooling state, 
but particularized for the fields n{z,t) and T{z,t). In this way, it follows that 

T{z-k,t-T)-T{z + X,t-T)p^T{z-X,t)-T{z + k,t) + T[T{z-k,t)C{z-X,t)-T{z + X,t)C{z + X,t)]. (43) 

Next, expansions around z, t are carried out as 

dT{z,t) 



T{z-X,t)-T{z + X,t)^-2X- 



dz 



T{z-X,t)i;{z-X,t)-T{z + X,t)i;{z + X,t)K-2X 
Since ^ oc nT^I'^ (see Eq. (6)), it is 



dT{z,t) 
dz 



■T{z,t) 



dz 



dz n dz 2T dz 



and substitution in Eq. (45) yields 

T{z-X,t)^{z-X,t)~T{z + X,t)^{z + X,t)f^-2X 
Use of Eqs. (44) and (47) into Eq. (43) leads to 

T{z-X,t-T)-T{z + X,t-T)^-2X 



3C(z,0 dT{z,t) ^ T{z,t)az, t) dn{z,t) 



dz 



dz 



1 + 



dT{z,t) 2XTT{z,t%{z,t) dn{z,t) 



dz n{z,t) dz 

Finally, when the above expression is introduced into Eq. (40), it gives the generalization of the Fourier law. 



1z 



dT dn 
dz ^ dz^ 



(44) 
(45) 

(46) 
(47) 

(48) 

(49) 



with the following expressions for the thermal heat conductivity and the diffusive heat conductivity transport coeffi- 
cients: 



3nvX 



3tC 



u = -vXtTC. 
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Define 



3nvX 



(50) 
(51) 

(52) 



This is the elastic limit of Eq. (50) if the (small) dependence of v on the inelasticity is neglected. Then, Eqs. (50) and 
(51) become 



where 



,c(T) = Ko(r)(i + ^ 

^l{T) = ^^^{T)T^, 



(53) 
(54) 

(55) 



has been introduced. 
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FIGURE 3. Temperature and density profiles for a two-dimensional system with a = 0.9 and ; 
is scaled with some arbitrary temperature T{A) and the density is measured in unit of a^. 



; No/W = 7. The temperature 



Therefore, the very elementary reasoning presented above leads to and expression for the heat flow in a dilute 
granular gas having the same structure as the one derived by using more exact and complex methods based on 
the Boltzmann equation or Unear response theories discussed in the previous section. To carry out a more detailed 
comparison with the explicit forms obtained by the latter methods, define a dimensionless cooUng rate by 

r^Y' ^^^^ 

In terms of it, Eqs. (53) and (54) read 

K{T) = KoiT)K*, n{T)=Ho{T)n*, (57) 

with 

K*^1+3C, M*=2r- (58) 

It is well known that mean free path reasonings leads to values of the elastic heat conductivity Kq that differs from those 
obtained by applying the Chapmann-Enskog procedure to the Boltzmann equation by a factor of the order of unity. 
This comparison will be not repeated here and attention will be restricted to the values of the reduced coefficients k*, 
and n*. 

The reduced cooling rate ^* has been defined in Eq. (56). Note that this is a property of the HCS and, therefore, 
is assumed to be known. Agreement between the definition of in Eq. (6) and ^* is obtained if T is chosen to be 
T = 2£/vQ{t). Moreover, keeping only up to order I — a, Eqs. (33) and (34) yield 

''*~l + j3T^o(«), M*-^^Co(«). (59) 
The agreement can be considered as satisfactory, taking into account the looseness of the arguments used. 



HEAT FLUX IN AN OPEN FLUIDIZED GRANULAR GAS 

When an open granular in presence of gravity is fluidized by means of a vibrating bottom plate, it reaches a steady 
state in which the energy supplied by the vibrating wall balances the energy dissipated in collisions. Assuming that 
there are gradients only in the z-direction, the Navier-Stokes equations of the systems are [25]: 



dp 

3- = -nmg, 



(60) 
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FIGURE 4. Reduced diffusive heat conductivity ^* for a dilute gas of inelastic hard disks (left side) and spheres (right side) as a 
function of the coefficient of normal restitution a. the symbols are simulation results and the soHd lines the theoretical predictions 
discussed in Sec. . the error bars are the mean square deviations of the values obtained in different runs, varying also the number of 
particles in the system 



^^+n=o, (61) 

where p = nT is the hydrodynamic pressure, ^ is the local cooling rate of the gas, and g is the intensity of the 
gravitational field in the direction of decreasing z- In ref. [25] it was shown that the temperature profile soltuion of the 
above equation with the corresponding boundary conditions has a minimum at a certain height Zm, increasing from 
there on. Therefore, beyond the temperature inversion the heat current flows from cold to hot, in a counterintuitive 
way. On the other hand, the density profile exhibits a maximum (at a different position). Using Eq. (60), the heat flux 
given in Eq. (22) can be rewritten as 



Qziz) = n£vo{t) 



(62) 



At the temperature minimum, this equation becomes 

qz{zm) = 2nmg£vo{t)ii. (63) 

Note that the existence of the minimum and, therefore, of the increase of the temperature with the height, is directly 
related with existence of the new term in the Fourier law. Eq. (63) has been used in numerical simulations to measure 
the exponent n and check the theoretical prediction derived from the Boltzmann equation [26]. The Direct Simulation 
Monte Carlo Methods was used in order to explote the symmetry of the system, rendering the numerical methods 
much more efficient [27]. Both two-dimensional and three dimensional gases where studied. As an example, the 
steady temperature and density profiles for a system of hard disks (d = 2), with a = 0.9, and No/W = 7 , where is 
the number of particles and W is the width of the system, are shown in Fig. 3. Other details can be found in ref. [26]. 
The lines are guides for the eye. Both the minimum in the temperature profile and the maximum in the density one are 
clearly identified. 



In Fig. 4 the coefficient jj.* = jxri/ KqT = iJ.2£nv(){T) / Kq, where K"o is the elastic limit of the thermal heat conductivity, 
is plotted as a function of the restitution coefficient 0£ for d = 2 and d = 3. The symbols are results from the simulations 
using Eq. (63) to compute Ji. For each value of a several values of A^^ have been considered. Besides, different 
realizations have been run in each case. The error bars in the figure are the standard deviation of the obtained values. 
The solid Unes are the theoretical prediction given by Eq. (34), which is indistinguishable from the results reported in 



refs, [7| and [22]. Equation (62) has also been used in ref. [28J to measure the transport coefficient ji, but restricting the 
analysis to the quasi-elastic region. Moreover, profiles very similar to those in Fig. 3 have been observed in experiments 
[10, 11, 12, 13], although in some of them it was not associated to the modified heat current discussed here. 

Finally, let stress that the purpose here was not to argue about the existence of the term proportional to the density 
gradient in the modified Fourier law, but to provide some insight into the physical mechanisms leading to it, and 
to illustrate some observable consequences of the new term. It is worth to emphasize the crucial role played by 
the inelasticity of colUsions and the associated energy dissipation in the modification of the Fourier law. It is the 
time dependence of the reference homogeneous distribution function rather than the presence of large temperature or 
density gradients that is the origin of the new term in the expression of the heat flux. 
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